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Abstract 

In this paper, we review the development of new 
reduced-order modeling techniques and discuss their 
applicability to various problems in computational 
physics. Emphasis is given to methods based on 
Volterra series representations and the proper orthog- 
onal decomposition. Results are reported for different 
nonlinear systems to provide dear examples of the 
construction and use of reduced- order models, partic- 
ularly in the multi-disciplinary field of computational 
aeroelastieity. Unsteady aerodynamic and aeroelastic 
behaviors of two-dimensional and three-dimensional 
geometries are described. Large increases in compu- 
tational efficiency are obtained through the use of 
reduced-order models, thereby justifying the initial 
computational expense of constructing these models 
and motivating their use for multi-disciplinary design 
analysis. 

Introduction 

Prior to the recent appearance of powerful digi- 
tal computers, it was necessary to construct models 
of physical behaviors that took advantage of existing 
analytical techniques or which involved numerical cal- 
culations with small numbers of degrees of freedom 
(DOFs). Now, partial differential equations, represen- 
tative of complex physics that were previously unob- 
tainable, can be discretized and integrated with nu- 
merical algorithms implemented on massive parallel 
supercomputers. Indeed, the simulation of nonlinear 
physical behaviors in even three space dimensions has 
become relatively commonplace; problems with mil- 
lions of DOFs can be routinely simulated, thereby 
allowing investigators to capture precisely important 
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phenomena. For example, as part of the Department 
of Energy’s Accelerated Strategic Computing Initia- 
tive (ASCI), the world’s fastest supercomputer (ASCI 
White - 12 trillion calculation per second) is being 
used in a shift from nuclear test-based methods to 
computation-based methods. 

In the absence of other tools and analysis, numerical 
simulation is often insufficient in the knowing represen- 
tation of complex physics. We see two main limitations 
of numerical simulation. First, while simulation can 
provide detailed time histories of discretized field vari- 
ables, such data may not readily imbue the investiga- 
tor with an increased level of understanding concern- 
ing the physics essential to a given phenomenon. As 
is true of physical experiment, careful analysis of the 
data is required to develop simpler models that can be 
used to predict important characteristics of system be- 
havior. Tiiis process can be impeded by the enormous 
size of computed datasets. Second, without the ded- 
ication of massive resources, numerical simulation of 
large systems remains far too computationally expen- 
sive to be used in various multi-disciplinary settings, 
including control model synthesis, multi-variable op- 
timization, and stability prediction. For example, in 
the field of computational fluid dynamics (CFD), codes 
for the simulation of turbulent, viscous flows in three 
space dimensions are often used to obtain j joint solu- 
tions, but less frequently used in related disciplines, 
such as aeroelastieity and shape optimization. Thus, 
there is a fundamental gap between the analysis fi- 
delity available to simulate an individual case and that 
practical for multi-discplinary analysis. 

Both limitations of numerical simulation suggest 
that computed data need to be distilled into lower or- 
der models that can serve as the basis for additional 
analysis. The intent in constructing such reduced- 
order models (ROMs) is twofold: to provide quantita- 
tively accurate descriptions of the dynamics of systems 
at a computational cost much lower than the original 
numerical model, and to provide a means by which 
system dynamics can be readily interpreted. We think 
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of a low-order model as a characterization of a phys- 
ical process, such that the essential behaviors of the 
process are captured with a relatively small number 
of DOFs. The use of point vortices to simulate the 
nonlinear dynamics of vorticity generating systems is 
a simple example of low-order modeling. ROMs are 
low-order models derived from some appropriate pro- 
jection of a full systems DOFs to a much smaller set 
that encapsulates most, if not all, of the systems funda- 
mental dynamics. Model accuracy typically depends 
on the number of retained DOFs and the convergence 
properties of the ROM. The reduction in computa- 
tional cost needed to solve the ROM is offset by a 
potential loss of accuracy and model robustness. 

The term “fidelity” is used to denote the degree to 
which a model captures the physics of a phenomenon 
of interest. Low fidelity implies that a computational 
model is missing key physical behaviors that render the 
model highly inaccurate in certain regimes, whereas 
high fidelity implies a broader range of model applica- 
bility. These terms alone are of ambiguous meaning, 
as they are dependent on the class of problems to 
which models are applied. For example, the linear 
potential equation can form a very accurate basis for 
computing loads in the subsonic regime, especially if 
corrected for viscous effects, but is a low-fidelity rep- 
resentation of flow behavior in the transonic regime, 
since leading-order (nonlinear) physical behavior is not 
properly modeled with this equation. Even in the sub- 
sonic regime, the linear potential equation may not be 
regarded to form the basis for a high-fidelity model if 
vehicles at large angles of attack are to be simulated. 
For a specified range of problem interest and physi- 
cal behaviors, high-fidelity models capture the relevant 
physical behaviors using validated techniques to ac- 
ceptable levels of accuracy. Thus, use of high-fidelity 
models usually leads to accurate solutions. However, 
application of a high-fidelity model is not sufficient for 
accuracy, owing to the need to execute properly the 
model on a computer (e.g., proper construction of a. 
grid and selection of a time step). 

One important trend driving the development of 
new ROM techniques is the increasing level of fidelity 
within multi-disciplinary analysis and design, which is 
a necessary consequence of the need for increased per- 
formance and reliability in many systems. The general 
purpose of reduced-order modeling is to lower the com- 
putational DOFs present in a numerical model while 
retaining the model’s fidelity, i.e., the model’s abil- 
ity to capture physics of interest. Point simulations 
using high-fidelity equation sets (e.g., Navier-Stokes 
equations) typically cannot be obtained fast enough 
to permit design. This situation will improve, but at 
the same time, it is likely that models of even greater 
fidelity will be desired as more complex interaction 
physics are accounted for in simulation. 

This paper reviews recent progress in the develop- 


ment, of reduced-order modeling techniques and their 
application to multi-disciplinary problems, particu- 
larly in the area of computational aeroolasticity. The 
scope of the review is limited to the Volterra theory 
of nonlinear systems and the proper orthogonal de- 
composition (POD), which have been applied to the 
high-fidelity analysis of aeroelastic configurations in 
two and three space dimensions, and which show great 
potential for continued practical use. The remainder 
of this Introduction is devoted to an overview of these 
methods, while later in the paper, we describe appli- 
cations of these methods beyond the field of aeroelas- 
t.icity. 

Volterra Theory of Nonlinear Systems 

Over several decades of aerodynamics research, sim- 
ple analytical models have given way to numerical 
descriptions of vehicle flight loads. This transition 
has been described by Silva , 1 which we summarize 
here. Early mathematical models of unsteady aero- 
dynamic response capitalized on the efficiency and 
power of superposition of scaled and shifted funda- 
mental responses, or convolution. Classical mod- 
els of two-dimensional airfoils in incompressible flow 2 
include Wagners function (response to a unit step 
variation in angle of attack), Kussncr’s function (re- 
sponse to a sharp-edged gust in incompressible flow), 
Theodorsen’s function (frequency response to sinu- 
soidal motion), and Sear’s function (frequency re- 
sponse to a sinusoidal gust). As geometric complexity 
increased, the analytical derivation of response func- 
tions was no longer practical and the numerical com- 
putation of linear unsteady aerodynamic responses in 
the frequency domain became the method of choice . 3 
And, when geometry- and/or flow-induced nonlineari- 
ties became significant in the aerodynamic response of 
configurations of recent interest, the nonlinear equa- 
tions were computed via time integration. 

The trend towards time-domain numerical analysis 
of the aerodynamic equations has revealed the dynam- 
ics of numerous important and complex phenomena, 
but has not provided a framework for the analysis of 
complex configurations without severe computational 
costs. Aeroelastic analyses involving coupling of the 
nonlinear aerodynamic equations with the linear struc- 
tural equations have been particularly costly to carry 
out. Post-processing of time transients at numer- 
ous flight conditions can be used to compute stability 
boundaries of the coupled system, but this approach 
has not been used extensively in industry due to the 
prohibitively high computational costs. 

Attempts to address the problem of high computa- 
tional cost include the development of transonic indi- 
cia! responses . 4,5,6 Transonic indicial (step) responses 
are responses due to a step excitation of a particu- 
lar input, such as angle of attack, about a transonic 
(or nonlinear) steady state condition. Neural ner- 
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works have also been used to develop nonlinear models 
of unsteady aerodynamics 7 and nonlinear models of 
maneuvers (using an experimental database ). 6 Neural 
networks and Vol terra series have some similarities, 
since each involves the characterization of a system 
via an input-output mapping . 8,7 In particular, there 
is a direct relationship between the weights of a neural 
network and the kernels of a Vo It err a series represen- 
tation for a particular system . 9 

A major difference between Volterra series and neu- 
ral networks is in the training effort. Neural networks 
can require a substantial training effort 6 while Volterra 
series neither require a training period or curve fitting 
for model construction. Also, Volterra kernels pro- 
vide a direct means for physical interpretation of a 
system’s response characteristics in the time and fre- 
quency domains. However, potential disadvantages of 
the Volterra theory method include input amplitude 
limitations related to convergence issues and the need 
for higher order kernels . 1 

Another approach to reducing the computational 
cost of aerodynamic analysis with CFD is to restrict 
attention to linearized dynamics. The response of 
the linearized system about a nonlinear steady-state 
condition can be obtained with a linear state-space 
representation of the system at that condition. In 
this form, the order of the state-space model can be 
reduced using various techniques . 10, 11 One method 
for building a linearized, low-order, frequency-domain 
model from CFD analysis is to apply the exponen- 
tial (Gaussian) pulse input . 12 This method is used to 
excite the aeroelastic system, one mode at a time, us- 
ing a broadband Gaussian pulse. The time-domain 
responses are transformed into the frequency domain 
to obtain a frequency-domain generalized aerodynamic 
force (GAF) influence coefficient matrix. These lin- 
earized GAFs can then be used in standard linear 
aeroelastic analyses . 13 R a veil, Levy and Karpel 14 ap- 
ply this method while replacing the Gaussian pulse 
with step and impulse excitations. The responses to 
these alternate excitations can then be transformed 
into a state-space form for direct use in other disci- 
plines such as controls or optimization . 15 We describe 
other frequency-domain representations in the context 
of the proper orthogonal decomposition. 

However, in order to develop robust and efficient, 
nonlinear CFD-based ROMs that are mathematically 
correct, a. rigorous method, well defined in the time 
and frequency domains for continuous- and discrete- 
time systems, is required. The Volterra theory of 
nonlinear systems fulfills these requirements. In par- 
ticular, tli is theory has found wide application in the 
field of nonlinear discrete- time systems 16 and nonlin- 
ear digital filters for telecommunications and image 
processing. 1 ' 

Application of nonlinear system theories, including 
the Volterra. theory, to modeling nonlinear unsteady 


aerodynamic responses has not been extensive. One 
approach modeling unsteady transonic aerodynamic 
responses is Ueda and Dowell’s 18 concept of describ- 
ing functions, which is a harmonic balance technqiue 
involving one harmonic. Tobak and Pearson 19 apply 
the continuous-time Volterra concept of functionals to 
indicial (step) aerodynamic responses to compute non- 
linear stability derivatives. Jenkins 20 also investigates 
the determination of nonlinear aerodynamic indicial 
responses and nonlinear stability derivatives using sim- 
ilar functional concepts. Stalford, Bauman, Garrett, 
and Herdman 21 develop Volterra models for simulating 
the behavior of a simplified nonlinear stall/post-stall 
aircraft model and the limit cycle oscillations of a sim- 
plified wing-rock model. In particular, they establish 
a straightforward analytical procedure for deriving the 
Volterra kernels from known nonlinear functions. 

A particular response from a CFD code may pro- 
vide information regarding the nonlinear aerodynamic 
behavior of a complex configuration due to a particu- 
lar input at a particular flight condition. It does not, 
however, provide general information regarding the be- 
havior of the configuration to a variation of the input, 
or the flight condition, or both. As a result, repeated 
use of the CFD code is necessary as input parameters 
and flight conditions are varied. A primary feature 
of the Volterra approach is the ability to characterize 
a linearized or nonlinear system using a small number 
of CFD-code analyses. Once characterized (via step or 
impulse responses of various orders), the functions can 
be implemented in a computationally efficient convo- 
lution scheme for prediction of responses to arbitrary 
inputs without the costly repeated use of the CFD 
code of interest. Characterization of the nonlinear re- 
sponse to an arbitrary input via the Volterra theory 
requires identification of the nonlinear Volterra kernels 
for a specified configuration and flight condition. 

The problem of Volterra kernel identification is 
addressed by many investigators, including Rugh , 22 
Clancy and Rugh , 23 Schetzen , 24 and more recently 
by Boyd, Tang, and Cliua 25 and Reiscnthel . 26 There 
are several ways of identifying Volterra kernels in 
the time and frequency domains that can be applied 
to continuous- or discrete-time systems. Tromp and 
Jenkins 2 ' use indicial (step) responses from a Navier- 
Stokes CFD code and a Laplace domain scheme to 
identify the first-order kernel of a pitch-oscillating air- 
foil. Rodriguez 28 generates realizations of state-affine 
systems, which are related to discrete-time Volterra 
kernels, for aeroelastic analyses. Assuming high- 
frequency response, Silva 29 introduces the concept, of 
discrete- time, aerodynamic impulse responses, or ker- 
nels, for a rectangular wing under linear (subsonic) 
and nonlinear (transonic) conditions. Silva 30 improves 
upon these results by extending the methodology to 
arbitrary input frequencies, resulting in the first iden- 
tification of discrete-time impulse responses of an aero- 


of 2G 


American Institute of Aeronautics and Astronautics Paper 2001 - 0853 



dynamic system. It should be noted that owing to 
separation of the input terms, Silva’s first approach 
had limited applicability for the identification of non- 
linear Volterra kernels, 30 a situation which has been 
resolved more recently. 1 

In his dissertation, Silva 1 discusses the then prevail- 
ing misconceptions concerning aerodynamic impulse 
responses, including the purported difficulty in com- 
puting such responses. These misconceptions primar- 
ily arise from fundamental differences between tradi- 
tional, continuous-time theories and modern discrete- 
time formulations. The appearance of discrete-time 
methods has great implications for aeroelasticity and 
aeroservoelasticity by providing a means for efficiently 
modeling nonlinear aerodynamics. In a similar fash- 
ion, other fields are likely to benefit from the coupling 
of large simulation codes and discrete-time response 
methods. 

The Proper Orthogonal Decomposition 

Before discussing the background of the proper or- 
thogonal decomposition, as applied to large, discrete 
systems, we will first summarize the role POD of- 
ten plays in computational physics. Reduced-order 
modeling with POD is essentially analysis by an em- 
pirical spectral method. With spectral methods, field 
variables are approximated using expansions involving 
chosen sets of basis functions. The governing equa- 
tions are manipulated to obtain sets of equations for 
the coefficients of the expansions that can be solved 
to predict the behavior of field variables in space and 
time. The POD is an alternative basis that is derived 
from a set, of system observations. In short, samples, or 
snapshots, of system behavior are used in a computa- 
tion of appropriate sets of basis functions to represent 
system variables. The POD is remarkable in that the 
selection of basis functions is not just appropriate, but 
optimal, in a sense to be described further in the Anal- 
ysis section. 

The need to obtain samples of system behavior to 
construct the POD-based ROM is both a strength 
and a weakness of the method. One strength is that 
models can be efficiently tuned to capture physics in 
a high-fidelity manner. Two noteworthy weaknesses 
are the need to compute samples with a high-order, 
high-fidelity method, and the possible lack of model 
robustness to changes in parameters that govern sys- 
tem behavior. Generally, the payoff in applying POD 
is quite high when, following an initial investment of 
computation, a compact ROM can be constructed that 
can be used many times in, sav, a multidisciplinary 
environment and which is valid over a useful range of 
system states. 

The POD basis, otherwise known as the Karhunen- 
Loeve basis, 31,32 is not new, but dates to the 1940s for 
continuous systems. Use of the POD is also known as 
principal-component analysis in the statistical litera- 


ture. 33 In fluid mechanics, the POD was first applied 
to the study of turbulence and the analysis of t urbulent 
flow data. 34,35 Numerous studies since then have em- 
ployed POD to characterize the turbulent properties 
and dominant, or coherent, structures of wall-bounded 
flows and free shear layers using experimental data. 
This work, and much other activity related to the 
POD is thoroughly reviewed by Berkooz, Holmes and 
Lumley, 36 where references are given for early applica- 
tions of POD in the fields of image processing, signal 
analysis, data compression, chemical engineering, and 
oceanography. Other references can be found in the 
fields of civil engineering 37,38 and structural dynam- 
ics 39,40 

More recently, computational data has been used 
in the construction of POD bases. For example, 
Moin and Moser 41 used data from a numerically simu- 
lated channel flow to compute characteristic structures 
within the channel flowfield. Sirovich also put. forth 
the method of snapshots (or strobes) 42 to ease the 
computational burden of obtaining the K-L basis for a 
discrete system. With the method of snapshots, eige- 
nanalysis of an M x M matrix is carried out, where 
M is the number of snapshots, instead of an N x A r 
matrix, where N >> M is the number of data points 
in a snapshot. This process will be described in the 
Analysis section. 

With the successful interpretation of large compu- 
tational data sets using POD, the technique was ex- 
tended to the dynamical modeling of various systems, 
including fluid-dynamical systems. Through this ap- 
proach, for example, fluid-dynamical systems are first 
simulated with CFD techniques, samples are taken, a 
POD is constructed, and then a set of low-order equa- 
tions is formulated in the POD basis, typically with a 
Galerkin projection, to study the dynamics of the sys- 
tem. Deane, Kevrekidis, Karniadakis and Orszag give 
an early example of this process as applied to flows 
through grooved channels and about circular cylin- 
ders. 43 In their noteworthy work, they successfully 
apply POD-based ROMs to the prediction of limit, - 
cycle behavior in these systems. Other applications 
of POD to the dynamic modeling of nonlinear heat 
transfer and fluid dynamic problems can be found else- 
where. 44,40 ’ 40, 4 ' 

POD-based ROMs are now being developed for the 
analysis of aeroelastic systems. Romanowski authored 
the first paper, which appeared quite recently, docu- 
menting reduction of the aeroelastic equations using a 
K-L basis. 48 His time-domain procedure involved the 
construction of a low-order model for the linearized 
dynamics of an airfoil with structural coupling about, 
nonlinear static states computed with the Euler equa- 
tions. Subsequent to this work, frequency-domain 
procedures have appeared that more efficiently com- 
pute POD basis functions for linearized aeroelastic 
systems. 49,50 Work is also underway in the extension 
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of POD techniques to the analysis of problems in non- 
linear aeroelasticity. Beran, Huttsell, Buxton, Noll 
and Osswald 51 proposed and tested a computational 
framework for computing static and dynamic behav- 
iors of nonlinear systems with POD- based ROMs that 
Pettit and Beran 52,53 have extended to treat the dis- 
crete Euler equations. This work also forms the basis 
for an ongoing study of nonlinear panel response in the 
transonic regime. 54 Dowell, Thomas and Hall are also 
using POD techqniques to investigate the limit-cycle 
oscillation of an airfoil with a nonlinear structural cou- 
pling in the transonic regime. 55 


Analysis 

In this section, we review important aspects of the 
analytical foundation of Volterra theory and the POD, 
as well as the application of reduced-order modeling 
techniques based on these methods to aerodynamic 
and aeroelastic systems. As described above, the POD 
is being applied in many different scientific and engi- 
neering disciplines, including aeroelasticity. While ref- 
erences are drawn from numerous sources, this review 
of POD analysis is not intended to be comprehensive, 
but is focused on recent work with connectivity to Air 
Force research activities. 

We will restrict our attention to solution vectors in 
an A r -dimensional, real, Euclidean space; these vectors 
will be written in boldface, along with the matrices 
used in vector equations. Time-dependent vectors of 
spatially discretized field variables, refered to as full- 
system vectors, are written as w (t), where t is time. 

Volterra Theory 

We begin by reviewing key features of the Volterra 
theory, as applied to time-invariant, nonlinear, 
continuous- and discrete-time systems. The liter- 
ature on Volterra theory is rich, including several 
texts; 00, 57,22,58 we follow the presentations of Silva 1,59 
and Raveh, Levy and Karpel 14 to capture issues re- 
lated to aeroelastic analysis. Furthermore, this sec- 
tion will concentrate on time- domain Volterra for- 
mulations, consistent with the implied application to 
time-domain, computational aeroelasticity methods; 
the foundations and applications of frequency-domain 
Volterra theory can be found elsewhere, 07, 22,08,25 

While one goal of this work is to document the appli- 
cability of Volterra to discrete-time systems (e.g., sys- 
tems arising in CFD), we first consider time-invariant, 
nonlinear, continuous-time, systems. Of interest is the 
response of the system about an initial state w(0) = 
W 0 due to an arbitrary input u(t) (we take u as a real, 
scalar input, such as pitch angle of an airfoil) for t > 0. 
As applied to these systems, Volterra theory yields the 
response 


+ f I h 2 (t - r\,t - r 2 )u(Ti)u(T 2 )dridr2 + 

3 0 Jo 


[ h„(t - T X ,..,t- r„)u(r 1 )..u(r n )dr 1 ..dr„. 

0 

(1) 

The Volterra series in expression (1) contains three 
classes of terms. The first is the steady-state term sat- 
isfying the initial condition, ho = Wo- Next is the 
first response term, / 0 °°hi(t - r)u(r)dr, where h is 
known as the first-order kernel (or the linear unit im- 
pulse response). The identification of the kernel h(ri) 
is based on measuring the response of the system to a 
unit impulse (Dirac delta function) at T\ ~ 0. Equa- 
tion (1) requires the system to be time invariant, so 
that the system responds in an identical manner (but 
translated in time) to an impulse at any T\ >0. The 
first response term represents the convolution of the 
first-order kernel with the system inputs for times be- 
tween 0 and f, where by causality, inputs beyond time 
t are excluded. Lastly are the higher order terms in- 
volving the second-order kernel, h 2 , and the n Ul -order 
kernels, h n . These terms do not all vanish when the 
system is nonlinear. 09 For example, identification of 
the second-order kernel is based on measuring the two- 
dimensional response of the system following impulse 
inputs at two different times. More will be said about 
kernel identification shortly. 

The convergence of the Volterra series is dependent, 
on input magnitude and the degree of system non- 
linearity. Boyd 60 shows that the convergence of the 
Volterra series cannot be guaranteed when the maxi- 
mum value of the input exceeds a critical value, which 
is system dependent. Of course, the issue of conver- 
gence is important, since the Volterra series must be 
truncated for analysis of practical systems. Silva 1,59 
and Raveh ct, ah 14 consider a weakly nonlinear for- 
mulation, where it is assumed that the Volterra series 
can be accurately truncated beyond the second-order 
term: 

w (t) = h 0 + / h i(t - r)u{r)dr 
Jo 

+ f f h 2 (t - n,t - r2)u(Ti)u(T2)riTi(ir 2 . ( 2 ) 
Jo Jo 

The assumption of a weakly nonlinear system is consis- 
tent with the emergence of limit-cycle oscillation of a 
2-D aeroelastic system in transonic flow through a su- 
percritical Hopf bifurcation. 51 For linear systems, only 
the first-order kernel is non- trivial, and there are no 
limitations on input amplitude. 

Silva 1 derives the first- and second-order kernels, 
which are presented here in final form in terms of var- 
ious response functions: 



w (*) = h 0 4- 


hi (t — t)u(t)(1t 
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= “(w l (T],r 2 ) -w 0 (ri) - W 0 (r 2 )). (4) 

In (3), w 0 (ti) is the time response of the system to 
a unit impulse applied at time 0 and w 2 (ri) is the 
time response of the system to an impulse of twice 
unit magnitude at time 0. These response functions 
represent the memory of the system. If the system is 
linear, then w 2 = 2w 0 and hi = w 0 , which is why 
the first-order kernel is refered to as the linear unit 
impulse response. The identification of the second- 
order kernel is more demanding, since it is dependent 
on two parameters. Assuming r 2 > t\ in (4), w 0 (r 2 ) 
is the response of the system to an impulse at time r 2 . 

Time is discretized with a set of time steps of equiv- 
alent size. Time levels arc indexed from 0 (time 0) to 
n (time f), and the evaluation of w at time level n is 
denoted by w[n]. The convolution in discrete time is 

N 

w[n] — h 0 + ^ hi [n - k]u[k] 
k = o 

AT JV 

+ ^2 ^2 - fci, n - k2]u[ki]u[k2\. (5) 

*1 =0 fc 2 =0 

It should be noted that an important conceptual 
breakthrough in the development and application of 
the discrete-time Volterra theory as a ROM technique 
is understanding the fundamental difference between a 
continuous- time unit impulse response and a discrete- 
time unit impulse response. 1 ^ 9 The continuous-time 
unit impulse response is a highly abstract function 
which suffers from a difficult, if not impossible, prac- 
tical (i.e., numerical) application. The discrete-time 
unit impulse response (also known as a unit sample 
response), on the other hand, is specifically designed 
for discrete-time (i.e., numerical) applications. The 
proof of this and details regarding the very powerful 
unit sample response can be found in any modern text 
on digital signal processing. 61 

The identification of linearized and nonlinear 
Volterra kernels is an essential step in the develop- 
ment of ROMs based on Volterra theory, but it is not 
the final step. Ultimately, these functional kernels can 
be transformed into linearized and nonlinear (bilinear) 
state-space systems that can be easily implemented 
into other disciplines such as controls and optimiza- 
tion. 22, 1 Currently, research is underway to develop 
these models, and results should be available soon. 15 

In addition, some very interesting and fundamen- 
tal research in the area of frequency-domain Volterra 
theory 62 and experimental applications of Volterra 
methods as applied to nonlinear aeroelastic problems 63 
continues. 

Proper Orthogonal Decomposition (POD) 

POD is a linear method for establishing an optimal 
basis, or modal decomposition, of an ensemble of con- 
tinuous or discrete functions. Detailed derivations of 


the POD and its properties are available elsewhere 64,65 
and not repeated herein. In our discussion of POD, 
M basis vectors are used to represent deviations of 
w(t) from a base solution, W 0 . These are written as 
{e 1 , e 2 , ..., e M }, and are referred to by many names, 
including POD vectors, 50 empirical eigenfunctions 64 
or, simply, modes. 64,52 For the sake of brevity, we 
shall use the term “modes” to denote the POD basis 
vectors. The modes are orthonormal 

iT j f 1 if i = j 

e 6 | 0 otherwise, ^ 

and computed in a manner to be described shortly. 
The modal decomposition of w using M modes is given 
by 

M 

w « W 0 4- y: Wje 1 = W 0 4- <£w, (7) 

i=i 

where <1> is an N x M matrix containing the or- 
dered set of modes, = !e 1 ,e 2 ,...,e A/ ] and w is 
an M-dimensional vector of modal amplitudes, w ~ 
[uq , t&jvf ]. As a time- varying function, w is ap- 

proximated by W 0 + 4>w(t). 

As stated by Holmes et, al., 64 “Linearity is the source 
of the [POD] method’s strengths as well as its limita- 
tions ...” The method is linear owing to the inde- 
pendence of the modes from the modal amplitudes, 
thereby allowing for the straightforward construction 
of reduced -order equation sets from the full equation 
sets following mode computation. 

The POD modes are constructed by first computing 
samples, or snapshots, of system behavior (solutions 
at different instants in time for dynamic problems, 
or equilibrium solutions at different parameter values 
for static problems) and storing these samples in a 
snapshot matrix, S. For now, we assume that M snap- 
shots are collected and column-wise collocated into the 
N x M snapshot matrix: S = [w l , w 2 , w w ]. By 
assumption, the snapshot matrix represents a random 
vector class of signals associated with the system. The 
basis provided by the POD, known as the Karhunen- 
Loeve 31,32 (or K-L) basis, minimizes the error in ap- 
proximating a member of this class with fewer than M 
basis vectors. This property of optimal convergence 
associated with the K-L basis is established in many 
works. 64, 65,66,67,50 The K-L basis can be readily com- 
puted by relating the mode matrix to the snapshot 
matrix through a transformation matrix V , 3> = SV, 
maximizing the projection of the snapshot matrix onto 
the POD basis. This leads to the eigenproblem 

S r SV = VA (8) 

for eigenvectors V and eigenvalues A — diag(\i). The 
eigenvalues are non-negative, since S T S is symmetric 
and positive semi-definite. Provided that the eigenvec- 
tors V are scaled to be orthonormal, V T V = I (I is the 
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identity matrix), the transformation formula * = SV 
yields * 7 <l> = A. Multiplying each e* by y/Xi yields 
an orthonormal set of modes, * T * = I, as originally 
specified. It should be noted that for their frequency 
domain analysis, Hall et al. construct complex modes 
from S /y SV = VA, where S ;/ is the complex conju- 
gate of the transpose of S. 

In practice, fewer than M modes are retained to 
simulate system behavior. These are selected based 
on the size of the modal eigenvectors. Simply put, the 
K-L basis for a subspace of dimension M r < M is ob- 
tained by retaining the modes associated with the M r 
largest eigenvalues computed in (8). It should also be 
noted that K-L theory establishes the K-L basis to be 
the eigenvectors of SS 7 /M, where SS 7 is the snap- 
shot covariance matrix. Manipulation of (8) yields 
SS 7 * m *A, and from the singular value decom- 
position of S, M eigenvalues are equivalent to those 
in A, while the other N — M eigenvalues vanish (for 
finite-dimensional problems of dimension A r ). 68 Com- 
putation of V followed by evaluation of * = S V much 
more efficiently yields the POD modes than explicit 
analysis of the covariance matrix. In the following, we 
will consider the number of modes retained to be a 
variable denoted by M that is less than equal to the 
number of snapshots in S. 

The techniques described below provide different 
means for obtaining reduced-order equations sets gov- 
erning w(t) in the POD subspace. There are several 
methods for accomplishing the projection, including 
the Galerkin projection, “subspace” projection (for 
linear and nonlinear systems), and collocation. We 
will explore the former two approaches herein. 

Gove mi ny Equ a ti ons 

We place the nonlinear and spatially discretized, 
full-system equations in first-order form: 

c/w 

— = R(w;A), (9) 

where w is an array of variables associated with inte- 
rior evaluation points (e.g., cell centers) throughout 
a computational domain, A is a free parameter (or 
set of parameters), R is an array of nonlinear func- 
tions of the discrete variables, and t is time. For the 
discrete Euler equations in two space dimensions and 
conservative form, w is a collocation of density, two 
components of momentum, and total energy, involving 
N = 4N P variables, where N p is the number of interior 
evaluation points. We also define an array of variables 
associated with boundary grid points (so-called ghost 
points), w*,, that are referenced in tire evaluation of 
R but not explicitly carried as variables. The 
boundary conditions are specified in general form as 

B(w, w b ) ~ f (£), (10) 

where B can be nonlinear in w and w&, and f is a time- 
dependent array representative of an evolving bound- 


ary state. When f vanishes, equilibrium states, W, of 
(9) exist and satisfy R(W; A) = 0 and B(W. Wt) = 0. 

Linear (Frequency Domain) POD Formulation 

Rapid progress has been made in the application of 
POD to aerodynamic and aeroelastic equation sets us- 
ing a linear POD formulation. 50 The basic approach 
is to develop POD-based ROMs for the linearized dy- 
namics about equilibrium solutions of the fully nonlin- 
ear equations, (9). For a given value of A (e.g., Mach 
number), the nonlinear base solution is computed 
with CFD methods accelerated for steady-state con- 
ditions. The governing equations are then linearized 
for periodic disturbances of small amplitude, placed 
in frequency-domain form, and solved using similarly 
accelerated CFD methods. Solutions of the linearized 
equations, gathered for a range of different frequencies, 
serve as snapshots in the construction of a POD-based, 
linear ROM. Furthermore, a linearized ROM represen- 
tative of the aerodynamics, can be attached to a set 
of nonlinear structural dynamic equations to form a 
compact, nonlinear aeroelastic model. 69,50 

Following Hall et ah, 50 a small disturbance to the 
equilibrium state is introduced: 

w(t; A) = W(A) +q r 7 "", (11) 


where u is the disturbance frequency, q is a distur- 
bance of small amplitude, j is the imaginary number. 
The disturbance is a response to forced oscillation 
about the equilibrium state at the boundary (e.g., a 
response in fluid velocity to the time rat.e-of-eliange of 
angle-of-attack for a pitching airfoil): 

B(w ' Wi,) = = ^ b( — ( 12 ) 

where b represents the type of forcing applied to the 
system. Introducing w* = W& + leads to a 

small-disturbance boundary equation (assuming the 
invertability of ^— ^): 


OB 

O C 16 

0 w 6 


OB 


(13) 


This equation is coupled to the small-disturbance form 
of the governing equation, 

uR 

Jq+7^ — q*=jMi (14) 

ow b 


where J = ^ is the system Jacobian (an N x N real 
matrix). Combining both sets of equations yields 


Aq = juj b. 


(15) 


5R3B _l 0B h = — — _1 b 

~ <9w (, Ov//, dw + ^ ' Owi, Owo 

The solution of (15), q, is a function of tlie forcing 
frequency, u, and the character and amplitude of the 
forcing as expressed through b. 
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It should be noted that (15) is based solely on the 
spatial discretization leading to R and does not ac- 
count for the iterative character of numerical schemes 
by which q is computed. For example, Hall et al. 
obtain an equation with structure equivalent to that 
of (15) using an explicit, cell-centered, finite- volume 
Godunov method, but derive an expression for a Lax- 
Wendroff scheme that contains additional terms which 
are second order in cuP° 

The POD is constructed using M solutions q*(f — 
1 , ...,M) of (15), for different frequencies and forcing 
conditions (e.g., airfoil pitch and airfoil plunge), as 
column entries in the complex snapshot matrix, S. 50 
As has been reported, accuracy can be retained while 
keeping M quite small; M is typically between 10 and 
100. Approximately half of the solutions need not be 
computed from (15), since solutions q at negative fre- 
quencies are complex conjugates of those at positive 
frequencies. The complex mode matrix is com- 
puted by first solving the complex form of (8), and 
then forming the product SV. To predict the time- 
dependent response about the equilibrium state, q is 
approximated b}' <l>q (q is the array of reduced-order 
variables) and substituted into the small-disturbance 
equation (15): A<£q = juj b. Hall et al. 50 project this 
equation onto the space spanned by 4> to obtain a 
reduced-order set of equations: 

Aq = ju b (A = $ h A$, b = . (16) 

Returning to the original system, (15), we see the ad- 
vantage of the POD formulation. If several solutions 
of (15) are required for different forcing conditions, the 
matrix A may be £W-decomposod or may be analyzed 
for cigenmodes that will dominate in the predicted re- 
sponse. Such “up-front” computations reduce the per 
unit computational cost of solutions beyond the first, 
but become impractical when N becomes sufficiently 
large, since the computational effort grows at a much 
faster rate than the number of equations. For example, 
on square computational domains, the decomposition 
of A grows as N 2 . 

With POD, the up-front cost is far less than that 
just described. The empirical eigenvectors are com- 
puted once in no more than 0(NM 2 ) operations (the 
product S H S). These eigenvectors are used in the 
construction of A, which also requires 0(NM 2 ) opera- 
tions (assuming A is not explicitly formed, but implied 
through the computation of A3q as suggested by Hall 
et al. 50 ). In practice, M is sufficiently small that the 
work is dominated by computation of the snapshots, a 
process requiring O(NM) operations. Furthermore, 
for a specified level of accuracy, M does not typi- 
cally grow with A r ; i.e., beyond a nominal threshold, 
grid refinement does not better capture low-frequency, 
high-energy structures. Once (16) is formed through 
a set of transformations involving the empirical eigen- 
vectors, many different cases (O(N)) can be examined 


at a commensurate cost, owing to the smallness of A 
(an M x M matrix). Accuracy of the approach is 
high, provided that a sufficient number of modes are 
retained, and provided that the forcing conditions are 
within the scope of the sampling process. 

The POD approach is well-suited to multi- 
disciplinary analysis involving repeated interactions 
between equation sets. A POD-based ROM can be 
used to simplify a computationally demanding equa- 
tion set so that it can be efficiently integrated with 
a simpler equation set. For example, Hall et al. ap- 
ply their analysis to the study of an isolated airfoil 
in transonic flow with pitch (a) and plunge ( h ) struc- 
tural coupling. 50 Their development is now summa- 
rized. Following collocation of the structural variables 
into the array h = (/i,a) r , the structural dynamic 
equations are expressed as 


'-10' 

d 

' h ' 


" 0 

i ] 


h 


0 

0 M 

dt 

h j 

+ j 

K 

°J 


h 

— 

. F ( w ) . 


(17) 

where M is a 2 x 2 matrix containing the airfoil mass, 
static imbalance and moment of inertia (about the 
elastic axis), K is the 2x2 stiffness matrix, F is an 
array representing the integration of the discrete flow- 
field into an applied force and moment, and h = ^). 

Here, (17) is placed in small-disturbance form as- 
suming that aeroelastic equilibrium is achieved when 
h = 0. such that h = hoe JU,f and h = . The 

force and moment function is written as F = Cqe*", 
where q is the flowfield disturbance captured by the 
aerodynamic equations. C is a sparse 2 x N ma- 
trix that represents the discrete force and moment 
integration; it is a function of the reduced velocity 
(V = Uryc/idabyJJi, where is the freestream veloc- 
ity, is the pitch natural frequency, b is the airfoil 
semi-span, and // is the airfoil-fluid mass ratio). Us- 
ing a POD-based ROM, suitably trained for pitch 
and plunge oscillations, the disturbance q is replaced 
by a set of reduced-order variables as shown above 1 : 
q = <£q, In small-disturbance form, (17) becomes 


' 0 

I 


ho 

+ juj 

-I 0 


ho 


0 

K 

0 


.ho J 

o mJ 


.ho. 


C$q 


(18) 


Pitch and plunge variables are linked to the aerody- 
namic disturbance problem by defining a sparse A r x 2 
transformation matrix B such that b = Bh. Thus, 
(16) becomes Aq — = 0, leading to a cou- 

pled set of M + 4 equations: 



7 A 

0 

0 " 



q ' 



0 

0 

I 



ho 



[ -C<£ K 

0 



h 0 . 



0 

0 

1 


r q 

1 


0 

-I 

0 


h 0 


0 

0 

M 


ho 
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Equation (19) represents a complex, general eigen- 
value problem (for 6 — juj) of small size that can be 
used to compute the stability properties of the aeroe- 
lastic system with great efficiency. For a given airfoil 
configuration, the Butter speed can be bracketed by 
systematically varying reduced velocity as a parame- 
ter until the eigenvalue with largest real part changes 
sign. Tli is approach can be contrasted with the di- 
rect approach of Morton and Reran, 70 ' 71,72 which has 
been successfully applied to the prediction of airfoil 
flutter speeds in the transonic regime. Their direct 
method does not reduce the number of degrees of free- 
dom in tin* solution array, and amounts to an implicit 
analysis of the eigensystem /iq = Jq, which is ex- 
panded to include the structural equations. However, 
the method does search for a single, conjugate pair of 
eigenvectors that becomes neutrally stable, and with 
this restriction, allows flutter speeds to be predicted 
at. a computational rate comparable to that of solving 
the nonlinear equations for the static base solution. 
As will be further reported in the Results section, the 
POD-based approach extends nicely into three dimen- 
sions, 73 while much additional work in this direction 
is required for the direct approach. 

Nonlinear POD Formulation - Subspace Projection 

The linear POD formulation described above pro- 
vides a practical means for assessing the linear stability 
and dynamics of complex, acroelastic systems at a very 
small fraction of the computational cost of full-system 
analysis. Specific examples will be given in the Results 
Section. In cases where the dynamics are nonlinear, 
or in which a reduced-order model of the nonlinear, 
static behavior is desired, a different class of methods 
is required. Two methods are described for treating 
nonlinear systems with POD: a subspace projection 
technique 51,52 and the Galerkin projection technique. 
The former is described in this section, while the latter 
is treated in the following section. The following devel- 
opment assumes that computed fields are reasonably 
smooth; the issue of field discontinuities is partially 
addressed later. 

Equation (9) is projected onto the subspace of 
reduced-order variables through a weight. eel-residual 
approach. 74 The dynamic residual, 7v, is defined as 

dw 

TC = — - R(w;A), 

which, when forced to vanish after weighting by each 
of the M modes, yields 

<R r - R(W 0 + <£w; A) j = 0. (20) 

With = I, (20) takes a form equivalent to that 
previously applied by Pettit and Reran: 52 

^ = $ t R(W 0 + ^w;A). (21) 


Equilibrium solutions of (21) satisfy the equation 
R ee $ r R(W 0 + $w; A) = 0. (22) 


This system of M nonlinear equations can be effi- 
ciently solved with the chord method following com- 
putation of the Jacobian, J = |§. 5i With this ap- 
proach, the Jacobian is numerically computed about 
a specified state, w°, and then frozen in the itera- 
tive procedure J(w°) (w 1 ^ 1 - w") = -R, where the 
superscript index denotes iteration. Only O(M) eval- 
uations of R and 3>w are necessary to construct J. 
For the results presented below, chord iterates are con- 
tinued until ||R|| < 10~ 3 or a lack of convergence is 
demonstrated. 

Unsteady solutions, w (t), of (21) can be time in- 
tegrated using either explicit or implicit techniques. 
As described further below, results have recently been 
presented 75 for integration with the second-order- 
accurate Crank-Nicolson scheme: 76 

F(w” +1 ) = w n+1 -w n - y ($ t R(W 0 + #w”;A) 

+ <f> T R(Wo + #w” +1 ; A)) = 0, (23) 

where the superscript now denotes time level. At each 
time step, the nonlinear system (23) may again be 
solved with the chord method. For weakly nonlinear 
systems, the Jacobian I - A* J(w)/2 can be evaluated 
once at t = 0; with stronger nonlinearities the Jacobian 
can be periodically updated at additional computa- 
tional expense. 

Sensitivities of equilibrium solutions to a can be 
very efficiently predicted with POD-based ROMs. 
Three formulas are relevant: the relationship between 
full-system and ROM sensitivities, ^ the 

definition of ROM residuals, R — <F T R, and a condi- 
tion on the ROM solution path 

</R — Jdw + *77— da = 0, (24) 

aa 


dw _ <9R 

da da 


(25) 


From these formulas, it follows that sensitivities satisfy 

(26) 


dw - j r 0R 
da aa 


The ROM is most beneficial for sensitivity analysis 
when there are several different parameters on which 
the system depends. In the procedure above, J can 
be computed and decomposed once (as true for the 
steady-state analysis), and then repeatedly used in 
evaluating the sensitivity formula (26) for each vari- 
able. Following the decomposition of J, the primary 
computational expense in evaluating (26) is calculat- 
ing ^7, which is variable dependent, but efficiently 
obtained. 
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The point at which (21) loses stability to time- 
oscillatory disturbances is a Hopf bifurcation point. 
The stability exchange occurs when a complex pair of 
eigenvalues of J has vanishing real part, while other 
real components of the eigenvalues of J are nega- 
tive . 51 For small A/, this process can be inexpensively 
and indirectly tracked through examination of all the 
eigenvalues of J. We define 7 such that 7 (a) is the 
real part of the eigenvalue of J(a) with largest real 
part; critical points of stability exchange, a*, satisfy 
7 ( 0 *) = 0 . For practical problems, where M may not 
be very small, critical points can be directly found in 
a manner similar to that applied to aerolestic systems 
by Beran and Morton . 70 The direct approach involves 
solving 7 — 0 through Newton’s method: 


<?7 

da 


(A" +1 «) = - 7 (a"), 


(27) 


where the correction s\ l '+ l a is typically relaxed with 
the parameter cdhopF & l ' +l = a 1 ' + t^hopf 

Nonlinear POD Formulation - Galerkin Projection 

In comparison to the subspace projection method, 
a more compact and efficient ROM can be obtained 
through Galerkin projection. This approach is the 
most common technique used to obtain nonlinear 
ROMs through the proper orthogonal decomposition, 
including applications involving equations of fluid mo- 
tion / 15,46 

The projection procedure is derived from the origi- 
nal partial differential equations, written as 


<9w 

— R( w ; A) = 0, 


(28) 


where w is now interpreted as a continuous field vari- 
able. As with the subspace projection method, the 
residual of (28) is forced to vanish after integrating 
the residual over the flow domain f!, weighted by each 
POD mode e* (& — 1,2,..., M): 


L 


<9w 

~dt 


- R(w; A) df! = 0. 


(29) 


The modes are obtained by: carrying out a numerical 
sampling process to obtain S in discrete form; solving a 
eigenvalue problem similar to ( 8 ); computing = SV, 
and scaling the modes to be orthonormal: 


L 

J e l e J dQ = | 


s'svr/n - va, 


1 if i=j 

0 otherwise. 


(30) 


( 31 ) 


As an example, Park and Lee apply a Galerkin pro- 
jection to the reduction of the Navier-Stokes equations 
in non- conservative form for the two-dimensional flow 
of an incompressible fluid . 40 They obtain a set of M 


ordinary differential equations for w following substi- 
tution of w(£) = *w(i) into the continuum equations 

/ nr>\ f C 7. -i 


(29) (for k = 1 , M): 


(hh 

dt 


M 


M M 


= a k + + IE ? bf Q ki m , (32) 


/= 1 


1=1 m = 1 


where a k , P*/, and Qkim are elements of singly, dou- 
bly and triply indexed arrays computed by integrating 
over terms involving the M retained modes. For ex- 
ample, one term appearing in the expression for Qki. m 
is ffl e k e l which is evaluated numerically, since 

the modes are available in discrete form. 

The main difference between the Galerkin and sub- 
space projections is in the evaluation of spatial deriva- 
tives. With the Galerkin projection, modes are in- 
serted into terms involving continuum derivatives that 
are integrated once in weighted fashion prior to nu- 
merical solution of the ROM. Alternatively, with the 
subspace projection, modes are dynamically inserted 
into the full, spatially discretized equations, which are 
then projected onto the ROM subspace. There are 
advantages and disadvantages associated with each ap- 
proach. For the Galerkin projection, distillation of the 
original equations into a fully reduced-order system 
(e.g., (32)) yields a representation that can be very 
efficiently solved for steady and unsteady problems. 
However, the low computational cost of solving the 
ROM is accompanied by a high cost in the construction 
of the ROM. Returning to the Navier-Stokes equa- 
tions for incompressible flow, the cost of computing 
the elements Qkim is 0(NM 3 ), which is impracti- 
cal if M grows too large. Other potential drawbacks 
include: the need for ROM reconstruction if modes 
are changed; the lack of full-system data with which 
to assess model error; the need to have a problem- 
specific procedure by which ROMs are constructed, 
and the difficulty of expressing some parameters (e.g., 
those representing variations in local properties) in 
the ROM equations. The main disadvantage of the 
subspace projection method is evaluation of the full- 
system array R at each iteration, which increases the 
per iteration computational cost by at least O(N) and 
requires software linkage between the procedures that 
integrate the ROM and evaluate R. However, this 
method is free from the drawbacks associated with the 
Galerkin projection method. 

Nonlinear Dynamics with Harmonic Balance 

The last technique to be described in this review is 
the harmonic balance (HB) formulation of Hall et al .' 7 
Their approach was first applied to the computation 
of time-periodic, nonlinear, viscous flows in 2-D turbo- 
machinery cascades, and is now being used to analyze 
the aeroelastic behavior of airfoils in transonic flow / 8 
The goal of the HB method is the efficient computation 
of time-periodic solutions of large, nonlinear systems. 
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Current methods based on time integration, for exam- 
ple those applied to aeroelastic computations, 79,80 can 
be very demanding owing to the need for retention of 
time accuracy. The requirement for small, global time 
steps and accurate integration over numerous cycles 
increases computational cost over steady* state analy- 
sis, for which large, local time steps can be used. With 
the HB method, the governing equations can be recast 
in a steady-state form that, accounts for the under- 
lying time periodicity of the solution, and which can 
be solved with pseudo-time-integration using standard 
acceleration techniques (local time stepping and multi- 
grid acceleration). 

There are several techniques, including harmonic 
balance, the finite-difference method, the shooting 
method, and the Poincare map method, that have 
been applied to a large number of of relatively low- 
order problems involving nonlinear oscillations. These 
techniques are described in numerous texts and papers 
(and the references within). 81, 82,83,81 The method of 
multiple scales has also been applied to the simulation 
of limit-cycle oscillation for an airfoil/flap configura- 
tion ill transonic flow. 85 

However, the HB method developed by Hall et ah 
is designed to treat nonlinear, aerodynamic problems 
of practical size for which there are large shock mo- 
tions.' 7 To describe this technique, we choose w to be 
a discretized field variable with A r degrees of freedom 
satisfying w(t) = w(f + T), where T is the period of 
the oscillation. The solution is expanded in a Fourier 
series in time, truncated to 2N^b -f 1 terms: 

Nh b 

w (t) = w n e jwnt . (33) 

>i=~Nhb 

The term n = 0 corresponds to the mean flow. Ex- 
pansion (33) can be substituted into the governing 
equation (9) to obtain 2A T //n + 1 equations for the 
vector coefficients in the Fourier expansion, w n , by 
collecting terms of like frequency. Using the nomen- 
clature of Hall et ah, these equations are written as 
S(w) = R(w;A), where in the unsteady term, S rep- 
resents a collocation of the Fourier coefficients, 

S = jwNw, (34) 

R is a collection of nonlinear terms arising from the 
residual array, R, N is a diagonal matrix containing 
the harmonic indices (— A r //^, ..., 0, .... Njm), and w is 
the set of vector coefficients 

w = (w -jV " fl , ..., w°, ... , w Nf/B ) T . (35) 

The evaluation of R, what Hall et al. call 
the harmonic fluxes, is computationally expensive 
{0{N Nfj B )) for the Euler equations and not easily 
extended to turbulent, viscous flows. 77 They propose 
an alternative harmonic balance formulation that is 


both simpler and more efficient. First, w is defined at 
2 Nub 4- 1 instants in time, evenly distributed about 
the periodic orbit, and collected into a single array w* : 


(w(0), w 


2 A r HB + 1 



2NhbT \\ 
2Nhb + 1 / / 


Then, the discrete Fourier transform operator, E, 
which is an NNhb x NNjib blocked matrix, is used to 
relate w* and w according to w = Ew*. Substitution 
of this expression into the equation S = ju : Nw = R 
provides 

jLjE _1 NEw* = R% (3G) 

where R* is the residual array evaluated at the 
2Nhb 4- 1 temporal points. The right-hand side of 
(36) reduces to a simple form, since when the gov- 
erning equations are in strong-conservation form, the 
product E~ ! E = I can be formed through the linear 
derivative expressions in R. 

Finally, Hall et al. introduce a pseudo-time r n by 
which (36) is integrated towards “steady state:” 

A f 

+ jwE -1 NEw* = R*(w*;A). (37) 

UT 

This step is an important benefit of the HB method, 
because it allows the time-dependent solution to 
be computed with existing, accelerated steady-state 
methods that need far fewer iterations than time- 
accurate, time- integration methods. In this formula- 
tion. the “unsteady” term is replaced by ju;E -1 NEw* , 
requiring 0(NN] I ^) operations to evaluate (e.g., the 
multiplication of Ew*), an improvement over the HB 
method in fully spectral form. Furthermore, Hall et. 
al. report that the evaluation of the fluxes, RA which 
requires O(NNhb) operations using standard tech- 
niques, dominates the cost of the numerical scheme." 

We comment that the HB technique does not involve 
a reduction in the number of variables arising from 
spatial discretization, and does not provide a model 
that is a compact representation of the full system. 
However, this technique does yield an efficient and 
low-order representation of the temporal variations of 
complex systems experiencing cyclic behavior in time. 
Also, a form of this HB method that allows the period 
T to be explicitly treated as an unknown (i.e., for an 
autonomous system) is currently being developed by 
Hall and his colleagues. 


Results 

In this section, we present results obtained with the 
Vol terra theory of nonlinear systems and the proper 
orthogonal decomposition. Attention is primarily re- 
stricted to problems in unsteady aerodynamics and 
aeroelasticity, but results in other fields are drawn 
upon to show the wide applicability of the techniques. 
Volterra first-order kernals are used to simulate the 
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Fig. 1 Volterra kernels for CFD analysis of RAE airfoil. Left: first-order kernel and effect of ID 

amplitudes. Right: First five components of the second-order kernel for plunge. 


unsteady response of viscous flows about driven air- 
foils, and the aeroelastic response of airfoils and wings. 
Higher accuracy can be achieved with second-order 
kernals, and the relative improvement is quantified for 
selected configurations. We also report on Duke Uni- 
versity’s application of POD to the rapid prediction of 
flutter boundaries of airfoils and wings in the transonic 
regime 50 ’ 55 and Air Force Research Laboratory’s in- 
vestigation of reduced-order forms of the discrete Euler 
equations. 52,54 The section is concluded with a sum- 
mary of results obtained with the harmonic balance 
technique of Hall et al. 77 

Volterra Series Analysis 

First,- and second-order kernels for the Navier-Stokes 
solution (with Spalart-Allmaras turbulence model) of 
an RAE airfoil in plunge at a transonic Mach num- 
ber using the CFL3D code 86 are presented in Fig. 1. 
On the left are two sets of first-order kernels due to 
two different sets of excitation amplitudes. Recall that 
the first-order kernel is computed using (3) and is the 
result of two pulse responses, one at a particular ampli- 
tude and the second at double the first amplitude. One 
of the first-order kernels shown in Fig. 1 was computed 
using excitation plunge amplitudes of w = 0.01 and 
w = 0.02, where w is a percent of the chord of the air- 
foil. The other first-order kernel was computed using 
excitation plunge amplitudes of w = 0.1 and w = 0.2. 
It is clear that the two kernels are not linearly related, 
demonstrating how the first-order kernel can capture 
amplitude-dependent nonlinear effects. On the right of 
Fig. 1 are five components of the second-order kernel 
for this case. The second-order kernel is more com- 
plicated because it is a two-dimensional function of 
time. The important point to be made is that this 
kernel is readily generated and its relatively smaller 
values (compared with the first-order kernel) and its 
rapid convergence indicate a small (but not negligible) 
level of nonlinearity present. This information might 


be used to determine that the first-order kernel may 
be sufficient to capture the dominant nonlinear effects. 
This point is demonstrated in Fig. 2 

Fig. 2 is a comparison of plunge responses for three 
different plunge amplitudes for the same configura- 
tion. Specifically, a comparison is made between the 
full CFD solution due to a sinusoidal plunging motion 
and that obtained using the first-order kernel from 
Fig. 1 (due to the larger excitation amplitudes). As 
can be seen, the plunge response obtained using the 
Volterra first-order kernel compares identically with 
the response obtained from the full CFD solution for 
the two smaller amplitude responses. The compari- 
son for the largest amplitude response (i.e., nonlinear) 
is very good as well, with a slight but noticeable dif- 
ference between the two results. The nonlinearity of 
the large- amplitude plunge responses is confirmed by 
linearly scaling the smallest amplitude (i.e., linear) re- 
sponse which, as shown in Fig. 2, cannot capture the 
amplitude-dependent nonlinear plunge dynamics seen 
at the larger amplitude. The turnaround time (“wall- 
clock”) for the full CFD solution was on the order of a 
day whereas the Volterra first-order solution was com- 
puted on a workstation in 30 seconds using digital con- 
volution. The initial cost of computing the first-order 
kernel was trivial as well due to the rapid convergence 
of the pulse responses. In fact, since each pulse (unit 
and double amplitudes) goes to zero in less than 100 
time steps, the responses were generated using a de- 
bugging mode option available on the supercomputer 
system used. Using this option, computations requir- 
ing less than 300 time steps are executed immediately, 
intended for debugging purposes. As a result, each 
pulse was computed within five minutes, resulting in 
a first-order kernel that was computed in about ten 
minutes. Of course, once the kernel has been com- 
puted, it can be used to predict the response to an 
arbitrary input (steady, any and all frequencies, ran- 
dom) of arbitrary length via digital convolution on a 
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workstation. Using this method, there is no need for 
the repeated, and costly, execution of the CFD code 
for different inputs. 

Normal Force 



Fig. 2 Comparison of actual nonlinear and first- 
order Volterra responses for three different plunge 
motions and a linear response for the largest mo- 
tion. 

Raveh. Levy and Karpel have recently applied 
the Volterra-based ROM approach to analysis of the 
AGARD 445. G wing. 14 They simulate the flow field 
around the wing using the EZNSS Euler/Navier- 
Stokes code. 8 ' This code provides a choice between 
two implicit algorithms, t he Beam- Wanning algorithm 
and the partially flux-vector splitting algorit hm of Ste- 
ger et al. Grid generation and inter-grid connectivity 
are handled using the Chimera approach. The code 
was enhanced with an elastic capability to compute 
trimmed maneuvers of elastic aircraft. 87 For the CFD 
computations, the flow field around the wing was eval- 
uated on a C-H type grid, with 193 points in the chord- 
wise direction along the wing and its wake, G5 grid 
points in the span wise direction, and 41 grid points 
along the normal direction. 

A process of mode- by- mode excitation, discussed 
previously, was performed for this wing using four 
elastic modes at a Mach number of 0.9G. The mode- 
bv-mode excitation technique provides the unsteady 
aerodynamic response in all four modes due to an ex- 
citation of one of the modes. In this fashion, a matrix 
of four- by- four functions (sixteen total) is developed. 
Two sets of excitation inputs were used: the discrete- 
time pulse input and the discrete-time step input. The 
cost of computing these functions is minimal due to the 
rapid convergence of these functions and it consists of 
only four code executions. Once these functions were 
defined, several full CFD solutions were generated that 
were due to various sinusoidal inputs at various fre- 
quencies. Shown in Fig. 3 is just one of these results 
for a 5 Hz input frequency, comparing the result ob- 
tained from the full CFD solution to that obtained 
via convolution of the step or pulse responses with 
a 5 Hz sinusoid. As can be seen, the comparison is 


exact to plotting accuracy for most of the responses. 
The full CFD solution, consisting of 8000 iterat ions re- 
quired approximately 24 hours on an SGI Origin 2000 
computer with 4 CPUs. By comparison, the Volterra- 
based ROM response shown required about a minute. 
Even including the upfront cost of computing the pulse 
(or step responses), the computational cost savings are 
significant. More importantly, the same pulse (or step) 
functions can now be used to predict the response of 
the aeroelastic system to any arbitrary input of any 
length. 

Shown in the left and middle portions of Fig. 4 is a 
comparison of linear and nonlinear GAFs for the first 
two modes of the AGARD 445. G Aeroelastic Wing at 
a Mach number of 0.9G. Nonlinear GAFs refers to 
the GAFs computed using the Volterra pulse-response 
technique about, a nonlinear steady-state value by ex- 
citing one mode at a time and obtaining the resultant 
responses in the other modes. The CFD results are 
compared with those using the ZAERO code for a 
purely linear case. Frequency-domain values were ob- 
tained by performing a convolution of several frequen- 
cies of interest with the computed CFD-based pulse 
responses. As can be seen, the comparison is reason- 
able and shows the small (but not negligible) nonlinear 
aerodynamic effects identified using the Volterra pulse- 
response technique. 

But rather than transforming the time-domain 
GAFs into the frequency domain, discrete-time, state- 
space systems can be created using the Volterra pulse 
responses directly. Presented in the right portion of 
Fig. 4 is a comparison of the pulse responses for the 
AGARD 445. G Aeroelastic Wing due to a unit pulse in 
the first mode. The CFD-based pulses (circles) com- 
pare exactly with the pulse responses obtained from 
a st ate- space system generated to model this system. 
The 32nd-order state-space system is a complete rep- 
resentation of the entire frequency spectrum of the 
unsteady aerodynamics defined by the GAF influence 
functions for the four aeroelastic modes of this wing. 
The pulse responses due to unit pulses in the second, 
third, and fourth mode are just as good as those shown 
in Fig. 4, but are not presented here for brevity. 

POD Analysis 

The proper orthogonal decomposition has been ap- 
plied to a variety of multidisciplinary problems in- 
volving the aerodynamic equations. In the remainder 
of the Results section, we summarize recent, findings 
arising from the reduction of system order through 
POD-based modeling. We first review application of 
the frequency-domain POD to aeroelastic systems in 
two and three space dimensions, and then describe 
progress in the analysis of nonlinear equation sets with 
POD, including those that exhibit limit-cycle oscilla- 
tion. This section closes with a summary of stability 
results obtained for a front-stage rotor in viscous flow 
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Fig. 3 Comparison of direct and convolved responses to sinusoidal excitation at 5 Hz (Mach 0.96). 


using the harmonic balance technique. 

Linear (Frequency Domain) POD Analysis 

Hall, Thomas and Dowell 50 report an application 
of their linear POD formulation to a two-dimensional 
aeroelastic configuration, and, in an ongoing study, 
Thomas, Dowell and Hall 73 extend this previous work 
to the analysis of three-dimensional wings. In 2-D, 
Hall et. al. develop POD-based ROMs for the NASA 
Ames Research Center NACA 64 AO 10 airfoil with a 
pitch and plunge structural model representative of 
a swept-wing section. 88 Base solutions of the Euler 
equations, Wo, are computed with a node-centered 
Lax-WendrofF scheme for freestream Mach numbers, 
Moo, between 0.7 and 1.0. A shock first becomes ev- 
ident in the base flow solution near the airfoil crest 
at about Mach 0.8. The aerodynamic equations are 
solved on grids of O topology; solution insensitivity 
to grid refinement is verified using a sequence of grids 
involving 65 x 33, 97 x 49 and 129 x 65 nodes. 

Reduced-order models of the flow equations are con- 
structed from solution snapshots resulting from two in- 
dependent airfoil movements about base states: pitch 
oscillation and plunge oscillation. Using the proce- 
dure detailed in the Analysis section, snapshots are 
computed from (15) for reduced frequencies (nondi- 
mensionalized by freestream velocity and airfoil chord) 
evenly distributed between -1 and 1. Following collec- 
tion of 41 snapshots. Hall et al. assemble ROMs with 


up to 41 retained modes. 50 A ROM is computed for 
each different Mach number (and base flow) examined. 
The aerodynamic equations in reduced-order form are 
then coupled with the structural equations, leading to 
the low-order, aeroelastic eigenproblem (19). 

Hall et al. 50 carry out very efficient analyses of 
(19) to construct flutter boundaries for the thickened 
NACA 64 A0 10 airfoil. Their results are summarized 
in Fig. 5. Flutter speeds predicted with POD-based 
ROMs compare well with those previously reported 
from transonic small-disturbance analyses, 88,89,90 and 
are used to document precisely the fold in the flut- 
ter boundary characteristic of the NACA 64 A0 10 air- 
foil. As observed by Hall et ah, the Mach number at 
which the fold occurs is underpredicted by the small- 
disturbance analyses relative to the POD-based Euler 
analysis. 50 Fig. 5 also illustrates that the POD results 
are well converged in the number of grid points used 
in the CFD computations and the number of modes 
retained in the ROM. For only a small set of Mach 
numbers do the results obtained with either the coars- 
est grid or the fewest number of retained modes deviate 
from the remaining results. In particular, the highly 
nonlinear behavior around Mach 0.9 is well defined in 
most cases. 

In work soon to be published, 73 Thomas et al. have 
extended the linear POD formulation to the analysis 
of the (weakened) AGARD 445.6 wing. Their in- 
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Fig. 4 Volterra results for the AGARD 445.6 wing at Mach 0.96. Left and Middle: Comparison of linear 
and nonlinear GAFs for the first two wing modes. Right: Comparison of pulse responses. 



Mach Number, M Mach Number, M Mach Number. M 

Fig. 5 Flutter speed of NASA Ames Research Center NACA 64A010 airfoil for various Mach numbers 
computed with frequency domain, POD-based ROMs. Left: results from 41-mode ROM (finest grid) 
versus other computational results; Middle: sensitivity of 41-mode ROM to grid refinement; Right: 

sensitivity of POD results (finest grid) to number of retained modes (from Hall et al. 50 with permission). 


vestigation addresses two potential concerns for the 
application of POD-based methods to practical aeroe- 
las tic configurations: the size of ROMs necessary to 
capture 3-D effects, and the sampling requirements as- 
sociated with a multiplicity of configuration-dependent 
structural modes. Thomas et al. compute base flow 
solutions of the Euler equations with a grid of 0-0 
topology containing 49 x 33 x 33 nodes, and construct 
flutter boundaries through POD analysis for Mach 
numbers between 0.5 and 1.141. They develop POD- 
based ROMs using two approaches. First, a 5 5- mode 
POD is built at each Mach number by computing solu- 
tions of (15) at reduced frequencies evenly distributed 
between 0 and 0.5 (conjugate solutions are associated 
with negative frequencies) for the first 5 natural modes 
of the wing structure. With this surprisingly small 
number of POD modes, Thomas et al. obtain results 
that are very consistent, in terms of flutter speed and 


frequency ratio, with those published by Lee- Rausch 
and Batina. 91 Thomas et al. also propose and demon- 
strate a promising technique for reducing the number 
of snapshots necessary to construct an effective POD- 
based ROM for an aeroelastic configuration. With this 
approach, only two snapshots are required for each 
natural mode of the structure, in addition to a set 
of “fundamental” modes, to construct the aeroelastic 
ROM. 

Nonlinear POD Analysis - Subspace Pi'ojection 

To study the application of POD to problems in 
aeroelasticity for which the reduced- order model is 
nonlinear, Pettit and Beran have first examined flow- 
field response to steady and unsteady changes in struc- 
tural shape. One such problem which has proved 
valuable is the response of inviscid How over a deform- 
ing panel in two space dimensions. 52,75 The flowfield 
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is situated above an infinite, panel that lies in the 
y = 0 coordinate plane, except for a segment between 
x = and .t — \ for which the panel shape is 
a parabolic arc defined by y s {z,t) = <$(*) (l -4x 2 ). 
For steady problems, the amplitude parameter 8 is 
a specified constant, while for unsteady problems, 
8{t) ~ Si sin(cj^) (1 - e“ at ) is assumed, where <$i, 
and a are constants. The length of the “bump” and 
the far-field velocity, , serve in the nondimension- 
alization of the aerodynamic analysis. A schematic of 
the bump and the 71 x 141 baseline grid are shown in 
Fig. 6. The flowfiold response is assumed to be gov- 

De formed 1 Vinci Segment 

- u (Shown M ore I >c lo rmetl Th an Actual ) 



Fig. 6 Top: Schematic of panel and coordinate 
system. Bottom: Baseline grid. 

erned by the Euler equations, which are discretized 
and solved with the upwind total variation diminishing 
scheme of Harten-Yee. 92,93 A transpiration condition 
is used by Pettit and Beran to satisfy approximately 
the boundary conditions on the bump surface while 
not requiring grid deformation. 

Pettit and Beran 53 conduct an analysis of flow 
changes in response to static changes in amplitude 
of the bump described above, and their preliminary 
results are summarized here. The algorithm for full- 
system analysis is validated by comparing computed 
results with those obtained using Cobalt 6 o, an un- 
structured, finite-volume algorithm for the Euler and 
Navier- Stokes equations 94 that has been validated ex- 
tensively. Differences in the predictions of the two 
techniques are small for variations of both Mach num- 
ber and bump amplitude in the range of interest. 

A single, reduced-order model for the steady-state 
bump is constructed from 26 full-system solutions 
computed at Mach 1.1, 1.15, 1.2, 1.25 and 1.3. At 
each Mach number, snapshots are computed at several 
different values of 8. For a given Mach number, there 
is a critical bump amplitude beyond which the shock 
attached to the leading edge of the bump detached, 
forming a bow shock. At sufficiently large values of 
8 and prior to shock detachment, response of the sys- 
tem is nonlinear to changes in 8. Once the bow shock 
forms, flow structure (i.e., shock position) becomes 
highly sensitive to additional changes in Mach num- 
ber and bump amplitude, a situation not well suited 


for POD analysis. Thus, most sampling and ROM 
application is limited to cases for which the flow is 
entirely supersonic. 

With 26 solutions, a total of 104 modes are com- 
puted. Pettit and Beran examine results when all 
modes are included in the analysis and when there 
is truncation to 60 and 40 modes, 53 Full-system so- 
lutions of (9) are time integrated to steady-state with 
2000-10000 function evaluations, depending on param- 
eter selection, and are initialized with uniform flow. 
While these cases were computed serially, they could 
also have been computed in parallel. Up to a point, 
sampling is a naturally parallel process; as a solution 
space is revealed, fill-in cases can be computed when 
resources become available. Equilibrium solutions of 
the POD-based ROM, satisfying (22), are computed 
along paths of constant Mach number and varying 8 , 
starting with the trivial solution at 8 = 0 and pro- 
gressing in increments of A<* = 0.001. At each point, 
ROM values are initialized using the previous solution, 
and solutions are typically computed in 3-10 function 
evaluations. Computing on the order of 100 solution 
points with the ROM typically requires fewer function 
evaluations than that necessary for a single solution of 
the full-system. 

Pettit and Beran interpreted the steady-state results 
using the minimum local Mach number, which typi- 
cally is observed at the bump leading edge. 53 Exami- 
nation of ot her flow variables did not alter their conclu- 
sions concerning the viability of the POD-based ROM. 
Solutions are compared in Fig. 7. Full-system solu- 
tions not used as ROM snapshots are also computed 
to obtain predictions of system behavior away from 
snapshot locations. Norms of the full- system residual 
are readily computed after each function evaluation 
in the subspace projection method and used to evalu- 
ate the quality of reduced-order solutions. Modeling of 
parameter-space subdomains where residual norms be- 
come unacceptably high can be improved by retaining 
different sets of modes or through model reconstruc- 
tion. In the figures just cited, only ROM solutions 
with a residual norm less than 0.6 are displayed. Pet- 
tit and Beran found this cutoff value to be consistent 
with good ROM solutions for the configuration being 
investigated. When 60 or all 104 modes are retained, 
ROM solutions are highly accurate and reproduce the 
nonlinear behavior evident in the response of the full 
system. Interestingly, the 60-mode ROM is somewhat 
superior to that of the full, or 104- mode, ROM. Pettit 
and Beran speculate that higher order modes in the 
104- mode ROM, which possess length scales on the 
order of the node spacing, can lightly pollute solutions 
with the present scheme, which makes no attempt to 
filter out such modes. However, in the steady-state 
analysis, modes are retained that are much smaller in 
magnitude (eigenvalues of S T S as small as 10“ 10 ) than 
that of the dynamic analysis. The fundamental diffi- 
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culty with retaining low-energy inodes in the dynamic 
analysis is not observed in the steady-state analysis. 
This is not a practical problem, since it is desirable 
to use fewer than the maximum number of modes. 
When the number of retained inodes is decreased to 40, 
the accuracy of the ROM is nominally degraded, and 
the number of computed points satisfying the residual 
norm cutoff is diminished. Most notably, solutions at 
Mach 1.1 do not meet, this criterion, which may not 
be surprising, since these points are on the boundary 
of the snapshot region and reflect the most nonlinear 
fluid dynamic behavior. 

Beran and Pettit have also computed unsteady solu- 
tions of the Euler equations for supersonic flow over an 
oscillating bump with an implicit, POD formulation. 
Their results were recently presented at the ICES ’2K 
Conference, 75 and are reported here to document the 
robustness and accuracy of the implicit technique de- 
scribed in the Analysis section. They consider a Mach 
1.2 flow over a bump with oscillation specified by the 
following baseline parameters: = 0.005, u = 1.0, 

and u — 3. The computational domain is discretized 
in a manner equivalent to that used in the steady- 
state analysis. A POD-based ROM is constructed for 
the baseline case by explicitly integrating the discrete 
equations, expressed as (9), and collecting snapshots. 
Integration is carried out for 20 time units (approxi- 
mately 4 cycles) with a time step of 0.01, the maximum 
value observed to permit stable integration. A to- 
tal of 200 snapshots, each representing a collocation 
of the conserved variables over the computational do- 
main, are collected at 10-iterat.e intervals. The initial 
state of the flowfield is specified to be uniform flow at 
freest ream conditions and is used in the definition of 

W 0 . 

Pettit and Derail 02 block the snapshot matrix so 
that the non- trivial elements of each column of S rep- 
resent only one of the conserved variables at a specified 
instant. In this manner, the 200 snapshots described 
above fill S, which is organized as a block diagonal 
matrix with a total of 800 columns. Each block is as- 
sociated with one of the conserved variables, as are 
each of the 800 computed modes. This approach is in- 
efficient for very large problems owing to the increased 
size of S, but leads to an adaptable framework for com- 
puting modes for each conserved variable. 

Following (8), ROMs are computed for between 6 
and 20 retained modes. With 20 inodes retained, the 
number of degrees of freedom is decreased by a fac- 
tor of 2000. Pettit and Derail observe that with as 
few as 8 modes retained, ROM integration yields very 
accurate results in comparison to full- system analy- 
sis for the case used to construct the reduced-order 
model. This is illustrated qualitatively in the left and 
middle portions of Fig. 8, where the structure (and 
magnitudes) of the density fields near the bump explic- 
itly computed with a 14-mode ROM and full-system 


analysis are nearly identical (shown at the end of the 
sampling period, t = 20). In the right portion of Fig. 8, 
time histories of pressure at the bump midpoint are re- 
ported for a 10-mode ROM implicitly computed with 
the Crank-Nicolson scheme. There it, is seen that im- 
plicit integration of the ROM accurately reconstructs 
the aerodynamic response, even using time steps 40 
times larger than needed for explicit integration of the 
full system. Results are shown for 5 subiterates in 
the Crank-Nicolson scheme; solutions can be obtained 
with 2 subiterates (without loss of accuracy) an order- 
of- magnitude faster than with full-system analysis. 

In a manner motivated by the frequency-domain 
analysis of Hall et ah, 50 Pettit and Derail, 53 construct 
hybrid ROMs combining snapshots computed for dif- 
ferent bump amplitudes and frequencies. These hybrid 
ROMs successfully reproduce aerodynamic responses 
for cases not explic.tly included in the sampling pro- 
cess. Pettit and Derail also report 52 ’ 53 for the bump 
problem that the subspace projection method yields 
numerically divergent results when the number of re- 
tained modes exceeds 19. Dy increasing the number 
of modes, the onset of divergence can be delayed, and 
the accuracy of the ROM solution increased prior to 
divergence. At present, the cause of this instability 
has not been identified; however, a violation of mass 
conservation in the system can be associated with the 
onset of divergence, which suggests a direction for fu- 
ture examination. 

Nonlinear POD Analysis of Limit-Cycle Oscillation 

To assess the applicability of POD-based ROMs to 
differential equations that exhibit limit-cycle oscilla- 
tion (LCO), Beran et al. 5i computed solutions of a 
tubular reactor, known to experience LCO, y5, 90 with 
the subspace projection technique. The governing 
equations are 

^ = Lw t - uhQIwo), (38) 

at 

A 

= Lw 2 - Pi (vl’2 - Pi) + PiWiQiw-t), (39) 

at 

Ls k&‘-h «(<«>« 

where Pe, , $ 2 , #$, T, and /x are specified param- 
eters. Equations (38) and (39) describe convection, 
diffusion and reaction within the reactor, and are re- 
ferred to as the CDR equations. The variables aq and 
W ‘2 represent concentration and temperature, respec- 
tively, and the parameter /x (the Danikohler number) 
determines the ability of the CDR equations to sus- 
tain LCO. The spatial domain is normalized; boundary 
conditions are applied at x = 0 and x = 1. Following 
spatial discretization of the equations and specifica- 
tion of suitable initial conditions, which are details 
described elsewhere, 51 the CDR equations take the 
form % =R(w;/i). 
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Fig. 7 Minimum local Mach number achieved by full-system and ROM solutions at selected freestream 
Mach numbers and bump amplitudes: Left: All 104 modes; Middle: 60 modes; Right: 40 modes. 



Fig. 8 Comparison of full-system and POD-based ROM predictions of aerodynamic response to bump 
oscillation (<b = 0.005). Left and Middle: Computed density contours at time £ = 20 (14-mode ROM). 
Right: Time histories of computed surface pressure at bump midpoint for different integrations of 10-mode 

ROM. 

The CDR system experiences a supercritical Hopf procedure described above, while LCO solutions are 

bifurcation at /C = (i = 0.16504, 51 which is accu- found with an explicit procedure like that applied to 

rately predicted by a POD-hased ROM. The stability the full system. 51 It should be noted that even with 

properties of the CDR system are shown in Fig. 9, the explicit procedure, dynamic solutions of the ROM 

where it is seen that stability of the equilibrium sys- can be obtained with time steps 50 times larger than 

tern is lost beyond the bifurcation point. Solutions 0.0005. This increase in allowable time step is a. con- 

are characterized by the maximum value of tempera- sequence of the absence of high-frequency, odd-even 

ture computed over the domain, T mAx . The ROM is modes in the POD-based ROM that would destabilize 

developed by sampling the CDR system as it evolves the numerical scheme. 51 

towards steady-state (0 < £ < 2.5) for a value of /i lead- 

ing to system stability: n° = 0.16 < //*. Following the Equilibrium solutions of the full system and the 
procedure described above, 8 modes are computed and R0M are observed to be in excellent agreement. As 

retained, representing a 15-fold reduction in problem s f n in Fi S- 9 agreement is nearly exact at 

size. Solutions of the full system are explicitly com- = 0.16, where the POD is constructed, and is ex- 

puted with time integration using a maximum time cellent for the remaining values of n shown. A slight 

step of 0.0005 (limited by numerical stability). Equi- inaccuracy is introduced at the Hopf point, where new 

librium solutions of the ROM are computed with the system behavior becomes available. Beyond the Hopf 

point, LCO amplitude is well predicted with the ROM. 
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Fig. 9 Comparison of full-system and ROM solutions for an 8-mode ROM trained at f.i° = 0.16. Left: 
Equilibrium and LCO branches. Middle: ROM system stability. Right: Sensitivity to // variation. 


Quantitative differences between ROM and full-sytem 
solutions can be reduced with nominal improvements 
in sampling (increasing the number of retained modes 
or increasing the sampling period ). 51 

The critical value of fi at which the CDR system 
loses stability is also very accurately predicted using 
reduced-order modeling. The variation of the stability 
parameter 7 computed with the ROM (cf. (27)) is 
shown in Fig. 9 (middle). Stability loss is observed 
at // = 0.16503, nearly the same value predicted with 
the full-system equations (// = 0.16504). The Hopf 
point is computed directly with Newton’s method in 11 
iterations (tJhopf = 0.5), starting with the equilibrium 
solution at (i — 0.16. 

Sensitivity of the variable T max to the bifurcation 
parameter is computed with the ROM formulation 
(26) at a = 0.164 using the 8 -mode ROM described 
above. As shown graphically and quantitatively in 
Fig. 9, the accuracy is excellent, with only a 1% dif- 
ference between the ROM and full-system results. 

Nonlinear POD Formulation with Shocks 

The application of POD to flows with shocks is 
a challenging endeavor, owing to the obvious diffi- 
culty of capturing movements of solution discontinu- 
ities with a fixed set of global modes. In preliminary 
work, the utility of ROM/POD for modeling the quasi- 
static movement of strong shock waves in a quasi-one- 
dimensional nozzle, assuming in viscid flow, is being 
assessed . 97 The location of the standing shock is var- 
ied by altering the the boundary conditions and/or the 
ratio of specific heats. It is observed that with straight- 
forward application of the nonlinear POD analysis de- 
scribed above, accurate modeling of shock movement 
requires an excessively large number of modes and data 
samples. Essentially, one snapshot is needed for every 
grid point location traversed by the moving shock. 

To improve the effectiveness of POD for problems 
of this type, a domain decomposition approach has 
been successfully developed as part of an ongoing in- 
vestigation . 97 In this approach, regions of the flow field 


that do not experience shocks are modeled with POD, 
while the flow-field region over which the shock moves, 
identified during the sampling procedure, is modeled 
with the full-system equations. See Fig. 10. With 
almost no degradation of accuracy, this methodology 
captures shock movement in the nozzle using an order- 
of-magnitude fewer degrees of freedom. An extension 
of these techniques is currently underway for a two- 
dimensional configuration where greater levels of order 
reduction are anticipated. 



Section I Section II Section III 



Fig. 10 Top: Nozzle with standing shock. Bottom: 
Decomposed domain for POD analysis. 

Nonlinear POD Analysis - Galerkin Projection 

Galerkin projection is being used by a number of in- 
vestigators to develop POD-based ROMs of low order 
for the study of various nonlinear phenomena. One 
noteworthy example is the work of Cizmas (Texas 
A&M University) and Palacios (San Diego State Uni- 
versity), which involves the construction of reduced 
order models for two-phase flow, heat transfer and 
combustion in dense or dilute fluid-solids flows. Their 
ongoing research project has several objectives: nu- 
merical generation of a database that includes spatio- 
temporal samples of system variables; modal decom po- 
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Fig. 11 Snapshots of gas velocity in the y-coordinate direction at 10 equally spaced time instants of a 
gas/solid mixture. 


sition of these variables through POD analysis; identi- 
fication and separation of dominant spatial structures 
during the evolution of the fluid and solids phases; 
Galerkin projection of the governing partial differen- 
tial equations onto POD-basis functions to produce a 
low-dimensional set of ordinary differential equations, 
and development of visualization tools for low-order 
models. 

Preliminary results obtained by Cizmas and Pala- 
cios demonstrate the ability of the POD to capture 
efficiently the energy content in a gas/solid mixture. 
The 2-D configuration is jet flow into a fluidized bed 
(40cm x 60cm) of 500 micron particles. The governing 
transport equations are much more complex than the 
Navier-Stokes equations; 3 gas and 3 solids equations 
comprise the set. Illustrated in Fig. 11 are ten snap- 
shots of the ^/-component of velocity taken at equal 
intervals in time that partially represent the ensemble 
of snapshots over all system variables. The snapshots 
are used in the computation of 20 POD modes for this 
component of velocity that account for 80% of the sys- 
tem energy (ratio of sum of retained POD eigenvalues 
to sum of all POD eigenvalues). Flow asymmetries 
about the vertical axis are captured in several re- 
tained modes. Their ongoing work, which is yielding 
promising results, is now focused on computing a low- 
dimensional set of ordinary differential equations that 
govern the fluid-solid flow. 

Transonic Aerodynamics with Harmonic Balance 

We now describe results obtained with the harmonic 
balance technique of Hall et al., 77 as summarized in the 
Analysis section. Hall et al. consider the flowfield near 
the tip of a front-stage rotor for a high-pressure com- 
pressor/ 7 The flow is modeled as two dimensional and 
the configuration is treated as an infinite blade row by 
imposing a suitable periodicity condition. The inlet 
Mach number and Reynolds number are specified to 
be 1.27 and the governing equations are the Navier- 
Stokes equations, closed with the Spalart-Allmaras 
turbulence model. 98 A baseline grid of H-O-H topology 
with 193x33 points is found to produce grid-converged 


results. The complexity of the steady-state flowfield 
is captured in Fig. 12, where shock-shock and shock- 
boundary layer interactions are observed to occur. 



Fig. 12 Computed Mach number contours for 
transonic viscous flow through front stage compres- 
sor rotor (from Hall et al. 77 with permission). Two 
instances of spatially periodic flowfield shown for 
clarity. 

Time-periodic unsteadiness is introduced by Hall et 
al. in the rotor flowfield bv oscillating the rotors in 
pitch about their mid chords with a reduced frequency 
of 1. While they report many results, we summa- 
rize their findings for the imaginary component of the 
first harmonic of the unsteady surface pressure, which 
has implications for aeroelastic stability. 77 Variations 
with respect to two parameters are examined: the 
pitch amplitude (a) and the interblade phase angle 
(a). Pressure is nondimensionalized by inlet dynamic 
pressure (q\ n \et), airfoil chord (c), and a (to better 
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display nonlinear effects). For a — 30 degrees, first 
harmonics of the unsteady pressure distributions are 
shown for different values of pitch amplitude and num- 
bers of retained harmonics (Nub) in Fig. 13. The 
first harmonic solution is rapidly convergent in A r ///?, 
with 5 harmonics being sufficient to capture all details 
of the unsteady pressure distribution. Nonlinearity 
arises with increasing a, and are particularly notice- 
able when a reaches 1 degree near a normalized airfoil 
surface position of 0.4. At this location, the unsteady 
pressure spike spreads out, a result of shock movement 
about this point. 77 




Fig. 13 First harmonic of unsteady pressure dis- 
tribution for front stage compressor rotor airfoils 
vibrating in pitch. Top: Variation of number of re- 
tained harmonics for pitch amplitude of one degree. 
Bottom: Variation of pitch amplitude for Nhb = 5 
(from Hall et al. 77 with permission). 

Variation of the imaginary part of the first harmonic 
of the pitching moment yields a more dramatic indi- 
cation of the role of nonlinearity. Blade stability is 


achieved when this component of the integrated pitch- 
ing moment is negative for all values of a.' 1 Results 
are displayed in Fig. 14, where it is seen that nonlin- 
ear aerodynamics serve to stabilize the flow when d is 
increased to about 1 degree. 



Fig. 14 Imaginary part of first harmonic of the 
blade pitching moment as a function of interblade 
phase angle (from Hall et al. 77 with permission). 

The HB scheme is reported to be highly efficient, 
an order of magnitude faster than conventional time- 
accurate, time-integration schemes. 77 Converged re- 
sults are found when 5 harmonics are retained, at 
a computational cost 7.5 times larger than that of 
steady-state analysis, per iteration. Convergence to a 
numerical solution is not obtained with the HB method 
when Nhb is increased to 7, a situation that is being 
further investigated at Duke University. 

Thomas, Dowell and Hall are also applying the HB 
scheme to the study of limit-cycle oscillation in air- 
foils with structural coupling in the transonic regime. 78 
Tli is work is important in unraveling the role of moving 
shocks in the phenomenon of limit-cycle oscillation. 

Concluding Remarks 

The development of ROMs based on the time- 
domain version of the Voiterra theory of nonlinear 
systems has been described, including continuous- and 
discrete- time versions of the theory. The basic ob- 
jective of the theory is the identification of linearized 
and nonlinear kernel functions that capture the dom- 
inant response features of a nonlinear system. This 
is, in fact, a nonlinear Green’s function method that 
provides a very natural and intuitive extension of well- 
understood linear phenomena into the nonlinear do- 
main. 

In the fields of unsteady aerodynamics and aeroe- 
lasticity, the use of influence coefficient functions, such 
as aerodynamic influence coefficient (AICs) functions 
and structural influence coefficient (SICs) functions 
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are routine. Since these influence coefficient functions 
are derived from Green’s function concepts, the exten- 
sion of these concepts into the nonlinear domain via 
the Volterra theory is quite natural and effective. In 
fact, these functions can be seen as a linear subset of 
a much broader nonlinear Vol terra functional space. 

While in its infancy for the analysis of large, spa- 
tially discrete systems, the Volterra method is cur- 
rently under continued development. One of the en- 
couraging applications presented is for a Navicr-Stokes 
solution of a 2-D RAE airfoil in transonic flow using 
the CFL3D code with the Spalart-Allmaras turbulence 
model. First- and second -order kernels were computed 
for plunging motions of the airfoil. It was shown 
that the first-order kernel (that captures amplitude- 
dependent nonlinearities) was used successfully to pre- 
dict the plunging response of the airfoil for large (non- 
linear) plunging motions at a minute fraction of the 
cost of the full CFD solution. Another important re- 
sult reviewed was an Euler solution of the AGARD 
445.6 Aeroelastic Wing recently computed by Raveh 
et al. 14 using the EZNSS CFD code. It was shown 
that with pulse or step responses the full linearized 
frequency-domain generalized aerodynamic force ma- 
trices were reproduced at a fraction of the time and 
cost than that required by standard techniques. 

The status of the Volterra-based ROM approach 
can be summarized as follows. The method has been 
used to show that discrete-time concepts, indeed digi- 
tal signal processing concepts such as unit pulses and 
step inputs, are directly applicable to CFD codes. 
The method has also been shown to be a higher- level 
generalization of the standard linear methods in use 
today. This is beneficial, because it means that indus- 
try experts do not need to restructure their analysis 
process in order to introduce Volterra-based methods 
into their design algorithms. In addition, the nature 
of the method is such that it requires minimal, if any, 
modification to the CFD code of interest. Most un- 
steady aerodynamic or aeroelastic CFD codes already 
have various excitation inputs (e.g., sinusoidal) and 
extension to a. Volterra-based ROM approach simply 
involves adding a pulse (or step) input to the suite 
of available inputs - the CFD code itself remains un- 
changed. 

As for the challenges associated with the Volterra- 
based ROM approach, there is much work to be done. 
An important issue that needs to be addressed is the 
issue of modal superposition with respect to nonlin- 
ear effects. Although it is clear that a mode-bv-mode 
excitation is a linearization of the aeroelastic process, 
it is important to understand the limitations of this 
approach. In addition, work continues on the develop- 
ment of a technique that provides simultaneous exci- 
tation to all modes, eliminating the linearization issue. 
Linearized state-space models are being developed us- 
ing the CFD-based pulse responses. These state-space 


models can be incorporated directly into control sys- 
tem analysis, for example. These state-space matrices 
also sidestep the need to transform time-domain CFD 
loads into the frequency-domain only to transform the 
frequency-domain loads back into the time domain via 
rational function approximations. Using the Volterra 
approach, time-domain CFD-based information goes 
directly into creating time-domain state-space matri- 
ces, a more efficient process. But the ultimate chal- 
lenge lies in the creation of nonlinear (bilinear) state- 
space matrices which are mathematically related to 
the Volterra kernels. Some work has been done in this 
area., but there is significantly more work that needs 
to be done. 

The results reviewed in this paper demonstrate the 
viability of POD -based reduced order models for rapid 
analysis of aerodynamic and aeroelastic problems. The 
frequency-domain approach of Hall et al. has previ- 
ously been shown 50 to predict accurately the linearized 
response of structurally coupled airfoils in transonic 
flow. As is typical in the use of POD-based ROMs, the 
computational cost of constructing low- order models 
through the frequency-domain approach is dominated 
by sampling (snapshot collection). Solutions are re- 
quired for discrete distributions of frequency in each 
of the char act erisitc modes of deformation (pitch and 
plunge in the case of the structurally supported air- 
foil). The cost of obtaining a large number of snap- 
shots can exceed that required to execute a small num- 
ber of general simulations. However, several important 
benefits in the application of POD to aeroelastic prob- 
lems are suggested by the work of Hall et al. for 2-D 
airfoils 50 and Thomas et al. for 3-D wings. 73 Following 
the one-time construction of a very low-order aeroe- 
lastic model, effects of parametric variations in the 
structural dynamics model can be rapidly assessed. 
Furthermore, the consequences of more fundamental 
changes to the structural model, such as freeplav,’ 5 
can be understood with greater clarity. The exten- 
sion to higher level multi-disciplinary applications, 
e.g., aeroservoelasticity or design with aeroelastic con- 
straints, may be made practical by the presence of 
low-order models. Lastly, the frequency domain POD 
represents an efficient and accurate compression of the 
salient, dynamic characteristics of the aeroelastic sys- 
tem. In summary, the efficiencies of POD-based ROMs 
are realized when system properties can be character- 
ized and observed once for ROM construction and then 
be allowed to vary in new wavs compatible with the 
previous observations. 

Application of the subspace projection method to 
the steady-state analysis of in viscid, compressible 
flow over a hump has yielded significant insight into 
the suitability of POD-based methods for nonlinear 
aerodynamic and aeroelastic analyses. By collecting 
a nominal number of solution samples over a two- 
parameter space, defined by bump amplitude and 
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freestream Mach number, a POD-based ROM was 
constructed that preserved the nonlinear transition of 
supersonic flow states towards transonic states. Not 
surprisingly, we observed that the range of validity (as 
measured by a residual norm of the full system) of 
the ROM decreased as the number of retained modes 
was decreased. Interestingly, validity was longest last- 
ing in the linear region of the parameter space (small 
amplitude for given Macli number), suggesting that 
low-energy modes were important in the accurate cap- 
ture of nonlinear behavior, as also seen in the dynamic 
analysis. Computation of ROM solutions required ap- 
proximately two orders of magnitude fewer function 
evaluations, which allowed a much more rapid and de- 
tailed exploration of the previously sampled solution 
space. We speculate that such ROMs are more reliable 
and useful than interpolation, although this should 
be documented. Aside from the potential advantages 
of data compression (i.e., keeping fewer modes than 
fully represented by the snapshot data), which could 
be significant in large systems characterized by large 
numbers of parameters, the ROM/POD provides a 
compact set of degrees of freedom that can be varied 
to evaluate sensitivities, optimal configurations, and 
system stability, in a manner based on the discrete 
equations of the full system. 

Unsteady aerodynamic solutions were also obtained 
with the subspace projection method at a significant 
computational savings over standard analysis. Re- 
sults were reported for the response of a supersonic 
flow over a bump, like that described in the steady 
analysis, whose amplitude varied sinusoidally in time. 
Following construction of a reduced-order model, the 
time-dependent character of the reduced-order system 
was accurately and efficiently computed with a sub- 
iterate form of the implicit Crank-Nicolson scheme. 
Time-dependent solutions of the ROM were computed 
an order of magnitude faster than full-system analysis. 
As true of the steady-state analysis, the computational 
cost associated with application of the POD-based 
ROM was dominated by the cost of data sampling used 
in ROM construction. The relative significance of the 
sampling cost can be minimized by constructing hy- 
brid ROMs that account for frequency and amplitude 
variations, and which are robust over a wider range of 
possible bump dynamics. 

Efficiency of the subspace projection technique was 
derived from two sources: decrease of the number of 
variables that characterize the system, and increase of 
allowable time step. In the bump problem, the number 
of variables was decreased by three orders of magni- 
tude (40,000 to about 10). The method was not three 
orders of magnitude faster, however, since evaluations 
of the full-system source term (R in (9)) were required 
with this approach. For steady analysis, full-system 
evaluations were employed in the construction of the 
ROM Jacobian, and in unsteady analysis, full-system 


evaluations were also necessary in the computation of 
dynamic residuals. With the POD-based ROM, com- 
putational work associated with implicit portions (i.e., 
left-hand sides) of system equations is virtually elimi- 
nated. Thus, the subspace projection method is partic- 
ularly well suited for implicit formulations of nonlinear 
problems, such as steady-state, sensitivity and bifur- 
cation analyses. For unsteady problems, it was also 
found that time steps allowed by POD-based ROMs 
were an order of magnitude larger than that allowed 
by explicit, full-system analysis. The current approach 
should be compared to the computation of a relevant 
viscous flow using a standard implicit technique to de- 
termine potential savings for a practical problem. 

Once sampling identifies the most energetic POD 
modes, other techniques are available with which the 
governing equations can be reduced in order; for ex- 
ample, Galerkin projection can be used to fully project 
the governing equations. The relative merits and de- 
merits of the Galerkin and subspace projection meth- 
ods are described in the Analysis section. Currently, 
the Galerkin approach is being used by Cizmas and 
Palacios to develop a small set of ordinary differential 
equations representative of jet flow in a solid/liquid 
mixture. A subspace projection method is also be- 
ing adapted by Lucia, King, Reran and Oxley 97 to 
treat a CFD problem for which the computational do- 
main is decomposed to isolate a moving shock. Other 
techniques, such as collocation, should be explored 
that may allow the POD modes to be used in a more 
efficient manner than subspace projection, but with 
perhaps greater flexibility than Galerkin projection. 

An alternative approach to POD based on har- 
monic balance has been proposed by Hall, Thomas 
and Clark 7 ' for the efficient computation of complex, 
time-periodic systems. With their technique, the re- 
sponse of a rotor flowfield to rotor pitch oscillation was 
accurately simulated, and behaviors related to shock 
movement and shock/boundary-layer interaction were 
captured. By using a low-order representation of sys- 
tem dynamics, the harmonic balance technique yielded 
response predictions approximately an order of mag- 
nitude faster than with traditional techniques. 

In closure, several methods have been described in 
this paper that offer new potential for the compu- 
tational analysis of large, nonlinear systems. These 
methods share a common reliance on existing numeri- 
cal techniques, and in this sense do not replace tradi- 
tional methods. Instead, reduced-order and harmonic- 
balance techniques provide existing methods with a 
higher level of algorithmic operation that enables 
more sophisticated computations. For example, a 
POD-based ROM of a discretized convection-diffusion- 
reaction (CDR) system was described and shown ca- 
pable of determining a variety of important charac- 
teristics of the nonlinear system, including nonlinear 
static behavior, bifurcation to limit-cycle behavior, 
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and sensitivity to changes in system parameters. The 
CDR problem serves as an analog for the study of 
the aeroelastic properties of a wing, including static 
analysis (e.g., determination of a control-surface re- 
versal speed), dynamic analysis (e.g.. prediction of a 
limit-cycle oscillation amplitude), bifurcation analysis 
(e.g., at what reduced velocities of the nonlinear sys- 
tem does flutter occur), and sensitivity analysis (e.g., 
how do aeroelastic behaviors depend on structural pa- 
rameters). Aeroelastic analysis of all these bahaviors 
in a manner that is useful for structural or aeroser- 
voelastic design is well beyond traditional methods. 
It is by answering more difficult questions, typically in 
the framework of multidisciplinary analysis, that ROM 
techniques become attractive, if not necessary. 

There are several challenges that need to be over- 
come before ROM methods can be routinely applied 
to practical problems. We group these difficulties into 
three categories: construction, generality, and accu- 
racy assessment. Which ROM method should be ap- 
plied to a particular problem will probably depend on 
the relative significance of these issues to the spec- 
ified problem. In the first category, work needs to 
he carried out to understand what response behav- 
iors should be included in the construction of ROMs 
to model robustly the response characteristics of sys- 
tems with large numbers of parameters. In other 
words, how much sampling is required for a partic- 
ular system? Generality of the ROM approach is 
also an important, issue. Is the approach readily or 
stubbornly extendable to different problems involving 
different simulation tools? Can the ROM approach 
function with modern, shock-capturing, CFD meth- 
ods that incorporate turbulence models and deforming 
meshes with unstructured or structured/overset con- 
nectivities? And finally, the accuracy of ROMs must 
be quantifiable for confident use. Systematic proce- 
dures must he developed (such as residual monitoring) 
to self-check ROM solutions and highlight conditions 
upon which ROMs should be re-constructed. 
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